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Abstract 



We investigate the bound states of the Yukawa potential V{r) = — Acxp(— ar)/r, us- 
ing different algorithms: solving the Schrodinger equation numerically and our Monte Carlo 
Hamiltonian approach. There is a critical a = ac, above which no bound state exists. We 
study the relation between ac and A for various angular momentum quantum number I, 
and find in atomic units, ac{l) = A[Ai exp(— //_Bi) + A2 exp{—l / B2)], with Ai = f.020(f8), 
Bi = 0.443(14), A2 = 0.170(17), and B2 = 2.490(180). 
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1 Introduction 

About 70 years ago, Yukawa proposed his famous meson theory to describe the interactions 
between nucleonsfl^, in which the potential is given by 



Here A and 1/a are respectively the strength and range of the nucleon force. Last few decades 
have seen a lot of ongoing analytical and numerical efforts [H 13 El IHl ^1 ^1 ^] on 
exploring the properties of the Yukawa potential. Up to very recently, this model still receives 
great attention |14[ I15[ I16| \17\ I18j . for it plays an important role not only in particle/nuclear 
physics, but also in many other branches: atomic physics, chemical physics, gravitational plasma 
physics, and solid-state physics. In plasma physics, it is known as the Debye-Hiickel potential; 
while in solid-state physics and atomic physics, it is named as the Thomas-Fermi or screened 
Coulomb potential. 
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Critical phenomena exist not only in Quantum Chromodynamics (QCD) at finite tempera- 
ture/chemical potential jl9l I2L)1 1211 122j . but also in the Yukawa potential in Quantum Mechanics 
(QM) [SlIlElIiniinillHISl- For a = 0, the Yukawa potential reduces to the Coulomb potential, 
and it is known to have infinite number of bound states. For a = oo, there is no interaction at 
all, and the system is completely free. Therefore, when a is non-zero, the Yukawa potential has 
some features very different from the Coulomb potential. One expects that the number of bound 
states is limited, because the interactions are screened. When a is sufficiently large, one expects 
that the bound states disappear. 

In this paper, we study systematically the bound states and the critical behavior of the Yukawa 
potential, using different algorithms: (1) numerical solution to differential equation |23j: (2) Monte 
Carlo Hamiltonian approach |24j. developed by some authors of the present article. The effects of 
the input parameters on the wave function and energy, and therefore on the critical behavior of 
the system are investigated in great details. We also propose a general relation among the critical 
a, A and the angular momentum I. 

The rest of paper is organized as follows. In Sec. [21 we describe briefly the numerical methods. 
In Sec. El we show the results. The critical behavior is discussed in Sec. 0]and the discussions 
are given in Sec. [S] 



2 Algorithms 

2.1 Schrodinger equation 

Because there is a spherical symmetry in the Yukawa potential Eq. , the Schrodinger equation 
is reduced to a radial one 



(Pu{r) 



+ 



u{r) = 0, 



where u{r) is the reduced wave function, and U{r) is a / dependent potential 



U{r) = V{r) + 



2mr^ 



(2) 



(3) 



Runge-Kutta and Numerov algorithms are two widely employed methods I2S1 for numerically 
solving such a Sturm-Liouville problem. Let us take the Numerov algorithm as an example. With 
q{r) = {E — U (r)) 2m/T)? , the Schrodinger equation 



d^u{r) 
dr"^ 



+ q{r)u{r) = 



is differentiated as 

where e is the integration step, and 



(4) 
(5) 



Cn+l 



1 + J^Qn+l, 



2 



5e^ 
b 

Cn-l = 1 + Y^Qn-l- 



(6) 



The boundary condition is 'u(O) = u{rmax) = 0. One can carry out the integration from r = and 
f = ^max, and match the wave function and the its derivative in the potential well region. Then we 
can combine the Sturm theorem and bisection method to find the eigenvalues and eigen-functions. 



2.2 Monte Carlo Hamiltonian 

Although most ID QM problems can well be studied by Runge-Kutta and Numerov algorithms, it 
is impossible to extend them to high dimensional QM or quantum field theory. Monte Carlo (MC) 
method with importance sampling 25 is able to compute high dimensional (and even "infinite" 
dimensional) path integrals. Unfortunately, using the standard Lagrangian MC technique, it is 
extremely difficult to estimate wave functions and spectrum of excited states. Let us take as an 
example a new type of hadrons predicted by QCD, the so-called glueballs j26ll?7| . which are bound 
states of gluons. Wave functions in conjunction with the energy spectrum contain more physical 
information on glueballs. 

Some years ago, we proposed a new approach |24j (named Monte Carlo Hamiltonian method or 
MCH) to investigate this problem. A lot of models [21 11111211001 ED El ESI EH m QM have been 
used to test the method. This method has also been applied to scalar field theories j35| l36 l l37 | l38j . 

We briefly mention the basic idea of MC Hamiltonian 24 ' . Starting from the Euclidean action 
Se = J dt (mr^/2 + U{r)), we consider the transition amplitude in imaginary time between time 
t = and t = T for all combinations of positions r^, £ {ri, . . . , r^v}. 



Mij{T) =< n\e-^^l\, >= J [dr]eM-SE[r]/n] 



(7) 



which is a path integral|3n] and can be computed by MC simulation with importance sampling |25j. 
at least in QM and the scalar models. The N x N transition matrix M can be diagonalized by a 
unitary transformation and used to construct the effective Hamiltonian j24j . Then we obtain the 
eigenvalues and eigen-functions of the effective Hamiltonian. 

Although MCH has been very successful applied to many local potentials, a special treatment 
has to be taken for QM with potentials having singularity at r = 0. The trick is to replace V{r) 
by 



Veffir) 



exp{2r/R) - 1 
exp(2r/i?) + f 



V{r). 



(8) 



Veff{r) is equal to V{r) in the R ^ limit. We first compute the spectrum and wave functions 
at some finite R, and then extrapolate the data to i? ^ 0. In Ref. 4^, we successfully applied 
this method to the hydrogen system. 
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Figure 1: Veff{r) as a function of r. 

3 Numerical results 

Botli the Numerov and Runge-Kutta algorithms can be directly applied to the Yukawa potential 
with very high precision. For MCH, as mentioned in Sec. [2 we have to introduce a cutoff R for 
Veff{r). Figure n shows the dependence of Veff{r) on R. For computing the transition matrix 
elements Eq. JT)), we choose the regular basis: 

rj=jxAr, j = l,2,...,N, (9) 

and let ^ = 1 , m = 1, Ar = 1.0, and N = 100. 

Figures I2I4I show the reduced wave function of the first three bound states at / = 0, A = 1 
and a = 0.1. These figures indicate that when — > the MCH data approach the Runge- 
Kutta/Numerov data. Table ^ lists some spectrum data for / = 0, A = 1 and a = 0.1, 0.2, 
0.3, 0.4. The error analysis follows the method described in Ref. [IHl- Within error bars, the 
MCH data agree with those obtained by the Runge-Kutta/Numerov algorithms. We have also 
confirmed that for sufficient large Vmax, which is the limit of integration, the results are stable 
and independent on Vmax- 



4 Critical Phenomena 

As mentioned in Sec. Q for a finite a, the number of bound states is finite. Let us first take the 
system with / = 0, A = 1 and a = 0.1 as an example. Figures El and El show the reduced wave 
function of the third excited state with r^nax — 

100 and 

f max — 200 respectively. The location of 
the lowest valley goes further and further away from the origin, with increasing limit of integration 
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Figure 2: Reduced wave function for the ground state at / = 0, and A = 1 and a = 0.1. 




Figure 3: Same as Fig. ^ but for the first excited state. 
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Figure 4: Same as Fig. [2 but for the second excited state. 



a 


Quantum number n 


Runge-Kutta/Numerov 


MCH {R 0) 


0.1 


1 


-0.4070 


-0.4051 ± 0.0612 




2 


-0.0499 


-0.0463 ± 0.0069 




3 


-0.0032 


-0.0031 ± 0.0008 


0.2 


1 


-0.3268 


-0.3345 ± 0.0604 




2 


-0.0121 


-0.0120 ± 0.0036 


0.3 


1 


-0.2576 


-0.2453 ± 0.0584 


0.4 


1 


-0.1984 


-0.2069 ± 0.0545 



Table 1: Comparison of some eigenvalues for Z = 0, A = 1, and some a, between the Rung 
Kutta/Numerov and MCH methods, where the MCH data have been extrapolated to the K = 
limit. 
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fmax- It will goes to r = oo when rmax oo. The eigenvalue E for the third excited state is 
8.1037 X 10-^ for r^ax = 100, 1.5533 x 10"^ for r^ax = 200, and 6.4491 x 10"^ for Vmax = 300. 
I.e., E goes to zero when rmax oo. Therefore the third excited state is no longer bounded. This 
also applies to higher states. 




20 40 60 80 100 



r 



Figure 5: Reduced wave function for the third excited state at / = 0, and A = 1 and a = 0.1, with 

1"max — 100. 

Table m shows the number of bound states for / = and various A and a. From the table, we 
find that when a is sufficiently large, the attractive potential is so screened so that the system 
no longer has a bound state. The critical value for a is named critical screening parameter ac- 
After fine scanning a and studying the dependence of u{r) on r^ax^ we obtain in Tab. Elthe value 
of ac for I = and various A. It is clear that ac satisfy the relation 

acil = 0) = 1.1906A. (10) 

This tells us that the physics is the same for all A, after rescaling the variables. 

We have obtained the result of Tab. Olby two criteria: (1) There is a unstable peak (or valley) 
in the ground wave function, with its location moving to oo; (2) The 

eigenvalue of the ground state becomes zero or positive as r^ax — > oo. 

The second criterion is easy to understand for the / = case. In such a case, U{r) = V{r) < 0. 
If the eigenvalue E > 0, then g(r) in Eq. (jlj becomes non-negative. According to QM, we know 
the state is unbounded. When I ^ 0, special care must be taken. Figures 1711^1 plot U{r) for / / 
and various a. When a « ac, as shown in Fig. [3 there is a well in U{r), so that some bound 
states with E < might exist; Of course, a state with ii^ > is unbounded and the particle can 
travel to arbitrary large distance, because q{r) is non-negative in the right hand side of the well. 
While a >> ac, as shown in Fig. IHl the effective interaction is repulsive, and no bound state 
could exist. 



7 




50 100 150 200 

r 



Figure 6: Same as Fig. [51 but with r^ax — 

200. 
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Table 2: Number of bound states at / = for various A and a. 



A 


1 


2 


3 


4 


5 


ac 


1.1906 


2.3812 


3.5718 


4.7624 


5.9530 



Table 3: Dependence of ac on A for / = 0. 
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U{r) 




Figure 9: U{r) for / 7^ when a » ac- 



When a is equal to or just above ac, at which U{r) behaves as Fig. |SJ the second term in 
Eq. (121) dominates at intermediate and larger r. Is the state with eigenvalue E = 0+ bounded? 
Without sufficiently large rmax, one would conclude that the state with E = 0+ to be bounded. 
However, with sufficiently large Vmax, Qii^) becomes non-negative at large r, and the particle can 
penetrate the well and goes to r — > 00 through tunnelling. This means the state with E = 0+ is 
actually not a bound state. 

One should look into the wave function in great detail. Let us take I = 4, A = 1, and 
a = 0.03135 cts cin Gxainple. U{t^ bchQ^ves like Fig. [HI When T^nax 

G [300,4267), the system 

seems to have a bounded ground state with E ~ 1.84 x 10~^. The reduced wave function of the 
ground state is shown in Fig. llOt The wave function and energy are stable within this range of 
Vfnax- However, when rmax ^ 4267, as shown in Fig. 1111 a peak in n(r) at large r appears. The 
location of this peak goes to r — s- cxd as increasing r^ax, at the same time, the energy of the state 
also tends to zero. Therefore, there is no ground state for ^ = 4, A = 1, and a = 0.03135. 
Table m lists the value of ac for various I, where ki is defined by 



acjl) 
A ■ 



Figure IT^ shows ki as a function of /. We can fit the data by 

ki = Ai exp 



-^) + ^2exp(--i-) 



(11) 



(12) 



with the fitting parameters given in Tab. |S1 For small /, the first term is dominant because 
Ai > Ai- While for larger the second term is dominant, because B\ < i?2- This implies that 
ac decays exponentially for large I. 
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Figure 10: Reduced wave function for the ground state with / = 4, A = 1, a = 0.03135 and 
rmax < 4267. 
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Figure 11: Same as Fig. E3 but with r^ax = 4267. 



11 
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4 


5 


k 


1.1906 


0.2202 


0.0913 


0.0498 


0.0313 


0.0215 


I 


6 


7 


8 


9 


10 




k 


0.0157 


0.0119 


0.0094 


0.0076 


0.0063 





Table 4: Dependence of ki on I. 




Figure 12: ki as a function of /. 



fitting coefficient 


Ai 


Bi 


A2 


B2 


data 


1.020 


0.443 


0.170 


2.490 


error 


0.018 


0.014 


0.017 


0.180 



Table 5: Fitting parameters in Eq. (|12|) . 
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5 Discussions 



In the preceding sections, we have systematically investigated the bound states and critical phe- 
nomena of the Yukawa potential, by means of different numerical algorithms. The results for the 
critical screening ac are consistent with the literature. We also obtain a new relation between 
ac and the angular momentum /. 

When judging whether a state is bounded, one must carefully study the effects of the boundary 
Tmax oil the wave function and energy. We believe such an experience might be very useful for 
numerical study of quantum theory, in particular at criticality. 

While the Runge-Kutta/Numerov algorithm works better in ID QM, it docs not work in 
quantum field theory. In this case, the Monte Carlo methods play an important role. The 
Monte Carlo Hamiltonian method has been tested as a useful method for computing the wave 
function and spectrum beyond the ground state, at least in QM and the scalar lattice field theory. 
Wc believe that the extension to QCD will allow to explore more useful information on strong 
interactions. 
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